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ABSTRACT 

We present a new modeling tool for plane- 
tary nebulae, based on 3D photoionization cal- 
culations. Our goal is to show that all the in- 
formation provided by observations, regarding 
kinematics and morphology, have to be consis- 
tently accounted for, in order to get a real in- 
sight of the object. Only 3D simulations offer 
this possibility. From models for two theoreti- 
cal PNe, we show that the enhancement in the 
equatorial zone observed in several PNe is not 
necessarily due to a density gradient, as usually 
interpreted. It is also shown that asymmetric ve- 
locity profiles often observed (e.g., Gesicki et al. 
1998) can be easily reproduced. Observations 
providing a better insight on the morphology of 
the PN are discussed. 



1. Introduction 

Planetary Nebulae (PNe) show different morphologies, 
but since the works of Kwok et al. (1978) and Balick 



(1987), it is possible to reproduce their observed shape ba- 



sically with two types of geometry: spherico-elliptical ones 
and bi-polar ones (also called butterfly), all of them axi- 
symmctrical in a first approximation. 

Observed morphological differences may result from (1) 
different orientation of the line-of-sight and of the axis of 
symmetry and (2) the age of the nebula. Morphological se- 
quences have been studied by Balick (1987), Pascoli (1990) 
and Zhang & Kwok (1998), among others. Evolutionary 
theories based on the generalized interacting stellar wind 
model (GISW, see Franck & Mellema 1994 and references 
therein) seem to confirm these interpretations. 

Information obtained from photometric imaging and in- 
termediate resolution spectra have been complemented by 
narrow band images (e.g., Sahai et al. 1997) and high res- 
olution spectra (e.g., Walton et al. 1990), providing infor- 
mation on the kinematics of the nebula. 

The various observations can not be correctly interpreted 
without detailed modeling, which must consistently repro- 
duce all observational data. For some years, most of the 
data have been analysed by GISW models (Soker & Livio 
1989; Mellema, Eulderink & Icke 1991; Icke, Balick & Frank 
1992), which account only for the gas dynamics. A more re- 
alistic model to analyse morphologies and kinematics of PNe 
was proposed by Frank & Mellema (1994). These 2D sim- 
ulations include radiation processes and ionization balance, 



while the gas cooling through line emission is obtained from 
polynomial approximations for the line emissivities (Balick 
et al. 1993). 

Here a new modeling tool for studying the morphology 
and kinematics of PNe is presented. Assuming a velocity 
field, the results of a 3D photoionization code (Gruenwald, 
Viegas & Broguiere 1997, hereinafter GVB) are used to gen- 
erate images to be directly confronted with various kinds of 
observations. We address general PNe issues, raised by the 
observations, as gas distribution and morphology, bright- 
ness enhancement at the equatorial ring, the asymmetric 
emission line profiles (e.g. Gesicki et al. 1998). 

The computational methods used to model planetary neb- 
ulae are presented in §2. Results for two theoretical objects 
are shown and discussed in §3. The conclusions appear in 
§4. 



2. Computational methods 
2.1. The 3D photoionization code 

The 3D photoionization code is described in GVB. The 
input parameters are the ionizing radiation spectrum, the 
gas chemical composition and the spatial density distribu- 
tion. The nebula is divided in a great number of cells. Inside 
each cell the physical conditions, which are obtained from 
the ionization equilibrium and thermal balance, are assumed 
homogeneous. For each cell, the output are the emissivity 
of atomic lines and the physical conditions of the gas (elec- 
tronic density, fractional abundances, electronic tempera- 
ture). The number of cells defines the numerical accuracy. 
It depends on the specific object under analysis, as well as 
on the computer memory. 

The effect of the diffuse radiation can be accounted for. 
However, in its present version, the required computational 
time is too high (few tens of the on-the-spot case). Thus, 
we assume here the on-the-spot approximation, using case 
B parameters of Pcquignot et al. (1991) for the hydro- 
gen recombination coefficient. The effect is small enough 
and do not affect the general conclusions presented here. 
The H/3 luminosity is underestimated by 30% or less when 
the on-the-spot approximation is assumed. The percentage 
depends on the stellar temperature: the higher the stellar 
temperature the bigger the underestimation of the H/3 lumi- 
nosity. The differences between the line intensities relative 
to H/3 is less than 20%, except for the weak lines (< 0.05 
H/3). 

Emission line images are obtained after a line-of-sight in- 
tegration. The calculations are made assuming a stationary 
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approximation. In order to analyse the kinematics, a ve- 
locity field is assumed after obtaining the photoionization 
results for each cell. The data cubes, containing the results 
generated by the 3D code for the cells, are read by an IDL 
(RSI) code which provides the tools for filtering, rotation, 
projection, statistics, velocity maps and echellograms, as 
described below. 

In this paper, only axi-symmetric nebulae will be stud- 
ied. In order to increase the number of cells used in the 
calculations (increasing the numerical accuracy), 1/8 of the 
nebula is modeled. For the calculated 1/8 of the nebula, the 
number of cells is 50 3 = 125 000. In this case, the ionizing 
source is placed at one of the corners of the large cube. In 
order to restored the whole nebula, with the ionizing source 
at the center, the remaining 7/8 are added up, using the 
calculated 1/8 after appropriate rotations. 

2.2. Rotations 

The nebula is defined in a XYZ coordinate system, where 
X is the axis of symmetry, the origin being at the ionizing 
source (center of the rebuilt cube). The rotations and pro- 
jections operate on a data cube of xyz coordinates. In the 
following, the line of sight used for the projections will be 
y. Before any rotation the origins, as well as the axes of the 
XYZ and xyz coordinate systems, are coincident. 

The rotations of the data cubes are performed in the xyz 
space, plane by plane, using the two-dimensional ROT func- 
tion of IDL, in a \/3 times greater cube to avoid losing in- 
formation. If needed, up to three rotations around the three 
axes of the cube can be performed, in any order. 

If rotations around more than one axis are necessary, it is 
easier if one starts with a rotation around an axis perpen- 
dicular to both the axis of symmetry of the nebula (X) and 
the axis along the line-of-sight (y). In order to compare the 
theoretical images to observations, a second rotation of the 
data cube may be needed. 

2.3. Emission line images 

For each emission line the 3D code provides the luminos- 
ity in each cell. The total luminosity (corresponding to the 
whole nebula), as well as the intensity corresponding to a 
given aperture, can be obtained and compared to observa- 
tional data. 

The projected image on the plane of the sky is obtained 
for a given emission line by integrating the emissivities of 
the cells lying along the line-of-sight which is taken per- 
perdicular to one of the cube faces. 



Maps of emission line intensity ratios can also be ob- 
tained, especially those used to obtain the PN physical con- 
ditions. These maps can be transformed into density or 
temperature maps using relationships (e.g., Alexander & 
Balick 1997) between the projected line intensity ratios and 
these quantities. 

2.4. Line profiles 

In order to obtain a line profile, a cube of velocities at the 
same resolution as the 3D code output is generated accord- 
ing to a predefined law. In this paper we assume a radial 
velocity field. Since the projection is made on the y-axis 
direction, only the corresponding component of the velocity 
is retained (V y (x, y, z)). For each emission line, a velocity 
profile is generated, for each line-of-sight, by 

, , x e x (x,y,z) -[ ^'.•.y) ]' 

^( W )=£^ (w) -e 

with: 

AV(v, x, y, z) = V y (x, y,z)-v 
£{x, y, z) = ] /v t 2 h (x,y,z) + VS 

V t h(x,y,z) = v / 2kT e (x,y,z)/Am H 

where e\(x,y, z) is the emissivity of a given line inside the 
xyz cell, Vt is the turbulent velocity and V t h(x, y, z) is the 
thermal velocity of the corresponding atom of mass num- 
ber A. The local electronic temperature T e (x,y, z) is con- 
sistently obtained from the 3D code. 

Velocity contour maps and P-V maps (echellograms) can 
be generated. After masking through an aperture and con- 
volving with an instrumental profile, the theoretical data 
can be compared to observations. 

2.5. Average physical quantities 

For each emission line, average physical quantities of the 
region where the line is produced can be obtained. Aver- 
age values for the electron temperature, density, hydrogen 
ionization, and velocity are calculated by weighting theses 
quantities by the line emissivity. The general assumed for- 
mula is: 

< Q >a= ^ — -, ^ 

where Q is replaced by T e , n e , H + / H or v. These quantities 
can be obtained for the whole nebula, as well as for lines-of- 
sight characterized by an xz position, summing on xyz or 
y, respectively. Notice that Peimbert (1967) defined some 
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of these quantities weighting by n e .n(X l ) instead of by the 
line emissivity. 

Dimensionless temperature fluctuations (t 2 ) and density 
fluctuations (n 2 ) can also be computed: 

^ _ T,(Q( x >v> z )- < Q >\) 2 -(\{x,y,z) 



Qt = 



Empirical methods for chemical abundance determination 
are usually based on these averaged quantities and fluc- 
tuations, defined for the whole nebula. They have been 
generally applied to point-to-point observations, assuming 
a spherical symmetry. In this case, however, the method is 
not always correct, even for spherically symmetric objects 
(Gruenwald & Viegas 1992). A real nebulae is not necessar- 
ily symmetric and homogeneous. Therefore, only 3D models 
can provide tools for improving the empirical methods. 

3. Application to planetary nebulae 
3.1. Models 

Two simple theoretical PNe will be assumed (A and B) 
with two possible geometries (1 and 2). Models A corre- 
spond to higher black-body temperature, higher density and 
higher luminosity than models B (Table Q). On the other 
hand, type 1 models have a spherical cavity surrounded by 
a shell with a density gradient; the density increases lin- 
early with the angle, from the pole to the equator. Type 2 
models have a prolate elliptical cavity surrounded by a con- 
stant density shell. These two geometries are two simple 
cases which can reproduce the bipolar brigthness enhance- 
ment usually shown by PN images. Surely, the geometry 
of real objects is determinated by hydrodynamics. Thus, a 
further improvement, out of the scope of this paper, would 
be to use the 3D photoionization code adopting a geometry 
obtained from a hydrodynamic code. 

The input parameters related to the density distribu- 
tion are given in Table The following notation is used: 
n# (cav) is the cavity density; n# (shell) is the shell density 
(type 2 models), or its range (type 1); Ri„„ er is the cavity 
radius (type 1) or the semi- axis of the elliptical cavity (type 
2). Models A are inspired on the canonical PN of the Paris 
workshop (Pequignot 1986, see also Ferland et al. 1995). 

The adopted chemical abundances are the same for the 
four models, and are taken from the canonical PN, i.e., 
H:1.0, Hc:0.1, C:3(-4), N:l(-4), 0:6(-4), Ne:1.5(-4), Mg:3(- 
5), Si:3(-5), S:1.5(-5), Cl:4(-7), Ar:6(-6), Fe:l(-7). 

For all models, the PN is radiation bound in all directions. 



3.2. Spectroscopical results 

Results for the 3D models Al, A2, Bl, and B2, corre- 
sponding to the whole nebula, are listed in Table Only 
line intensities of the main emission lines, relative to H/3, are 
given, as some average quantities obtained from line ratios, 
using Alexander & Balick (1997) fits or the formulae given 



in §2.5 (lower part of the table). 

The differences between the results of Al and A2 models 
are not significant. For both models, the results for line in- 
tensities are within the values obtained for the canonical PN 
(spherical symmetry and constant density) with the various 
ID photoionizing codes (Ferland et al. 1995). Notice that 
our results were obtained using the on-the-spot aproxima- 



tion (see §2.1). From the point of view of the emission-line 
spectra, as well as of the derived physical conditions, no 
conclusion about the nebula geometry can be drawn, since 
both geometries lead to similar spectroscopical results. The 
same is true for a planetary nebulae with a lower stellar 
temperature, lower luminosity and density, as seen from the 
similarity between the results of models Bl and B2. 

For both (A and B) types of nebula, the determination 
of temperature and density, from the line ratios or from the 
formulae (§2.5), are consistent within 10%. The electronic 



temperatures obtained from the [OIII] and [Nil] line ratios 
are different because these lines come from different regions. 

EDITOR: PLACE TABLE [l] HERE. 



3.3. Imaging 

Emission line images often show a brightness enhance- 
ment which is usually interpreted as the signature of a den- 



sity gradient in the equatorial belt of PNe (Weedman 1968 
Khan & West 1985; Aaquist & Kwok 1996; [Zhang fc Kwok 



19981 ). Thc images of H/3, [O III]A5007, He IIA4686, and [N 
II]A6583 obtained from models A and B are given, respec- 
tively, in Figs. |l| and ||. The nebulae are projected along 
the Y axis (no rotation) . The top four panels correspond to 
type 1 models, while the botton four ones to type 2 models. 
Notice that type 1 and 2 models have different geometries, 
which might reproduce the observed enhancement. In type 
1 models, the gas density increases linearly with the angle 
from the pole to the equator, and is distributed around an 
inner spherical cavity. On the other hand, in type 2 models 
the gas is uniformly distributed around a prolate elliptical 
cavity. As suggested by the figures, the models are able to 
generate the brightness enhancement. Notice that the line 
intensity enhancement is due to the higher density in the 
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equator for type 1 models, whereas, for type 2 models, it 
appears because the density of the ionizing photons reach- 
ing the inner edge of the gas is higher in the direction along 
the minor axis. 

However, the real geometry can not be easily inferred 
from the shape suggested by the images. In fact, model 
Al images show an ellipsoidal shape, while for model A2 
the shape suggests a spherical distribution, which is the op- 
posite of the adopted shapes for the inner cavities. This 
contrasting behavior is less pronounced in models B, with 
lower stellar temperature and luminosity, as well as with 
lower density, than models A. The image shape given by 
model A2 is defined by the ionizing radiation dilution, in- 
dependently of the gas distribution. On the other hand, in 
model B2 the main factor defining the image shape is the 
shape of the inner cavity. 

From these results, we see that narrow band imaging 
may not differentiate the nebular geometry as previously 
assumed. In the next section we discuss a possible geome- 
try indicator. 

EDITOR: PLACE FIGURE | HERE. 
EDITOR: PLACE FIGURE HERE. 



3.4. Density maps 

As shown above, spectroscopic results and emission line 
imaging do not provide enough contraints on the nebula 
density distribution. The density indicator generally used 
is the [SII] line intensity ratio. The density maps obtained 
from the [S II] ratio (see §2.3) are shown in Fig. || for models 
Bl (top panel) and B2 (bottom panel). As expected, the 
density gradient adopted for model Bl is clearly perceived. 

EDITOR: PLACE FIGURE H HERE. 



3.5. Velocity profiles 

High resolution observation of PNe emission lines often 
show double-peaked profiles, interpreted as the emission of 
a geometrically thin expanding shell. Usually, the two peaks 
are not symmetric (see Gesicki et al. 1996, 1998 for seven 



PNe, and [Dudziak et al. 1999| for NGC3132). If the asym- 
metry were due to local extinction, the blue peak should 
always be more intense than the red one. Since this is not 
the case, another explanation for the asymmetry must be 
found. 



In the following, the effect on the line profiles, due to 
geometry and orientation, will be discussed. We show that 
the asymmetry in the line profiles can be reproduced using a 
3D photoionization code and the numerical tools described 
above (§2). In order to illustrate these results, one case is 
presented, using model Al with a given velocity law. A lin- 
ear law, as used by Sahu & Desai (1986), is quite natural. In 
the following we adopted V — a x r/r , with a — 20 km/s 
and r a = 2.4xl0 17 cm (maximum of the [O III] emissivity). 
The velocity profiles obtained after the rotation of the neb- 
ula by 50° around the z-axis and projection along the y axis 
are shown in Fig. ^. Each profile corresponds to a synthetic 
observation at a different position in the nebula. The as- 
sumed positions and aperture are shown in the upper panel, 
over the [O III] image. A turbulent velocity of 2 km/s was 
taken into account. Such value for the turbulent velocity is 
small, so that the adopted velocity field dominates the line 
profile. No convolution through an instrumental profile was 
applied. 

As seen in the figure, the profile changes from one re- 
gion to another. Depending on the position, asymmetric 
profiles are obtained, including asymmetric double-peaked 
profiles. For positions at the right side of the nebula the 
double peaked profiles, when asymmetric, show the blue 
peak more intense than the red one. On the other hand, 
due to the rotation around the z-axis, a red asymmetry 
(not shown in the figure) would be obtained for positions 
on the left side of the image. 

The mean velocity (as defined in § |2.5|) of the [O III] zone 
is < V >5oo7 = 22.1 km/s, which is comparable to what 
could be found from the peak separation on the profile cor- 
responding to the synthetic observation of the central part 
of the nebula. 

EDITOR: PLACE FIGURE | HERE. 

Other velocity laws, suggested by hydrodynamic calcula- 
tions, could be adopted when discussing observations from 
real objects. 



3.6. PV diagrams 

More accurate than line profiles observed at discrete po- 
sitions, Position- Velocity (PV) diagrams (also called echel- 
lograms) provide an effective diagnostic for the geometry 
of PNe. In the case of line tilts (Guerrero et al. 1997), 
PV-diagrams give constraints on the inclination of the PN 



(Marston et al. 1998). The observations and the 2D model 



of NGC 650-1 (Bryce et al. 1996) is another example of the 
use of PV-diagrams to deduce the geometrical properties of 
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PNe. Such a model is not detailed enough to reproduce the 
differences between the observed emission of [N II]A6584 
and [O III]A5007. A 3D photoionization code is needed. 
PV-diagrams obtained from our model Al rotated by 50° 
around the z axis and projected on the y direction are shown 
m Fig. § The upper-left panel shows the [O 111]A5007 image 
of the nebula. The upper-right and lower-left panels show 
PV-diagrams obtained through vertical and horizontal cen- 
tered slits, respectively. The lower-right panel illustrates 
iso-velocity contours, for V = 9 km/s. The z-PV-diagram 
shows the signature of the tilt of the PN around the z-axis, 
while the absence of an asymmetry in the x-PV diagram in- 
dicates that the Z axis of the nebula is perpendicular to the 
line of sight and coincides with the z-axis. Such simulation 
is available from the 3D code for all the emission lines and is 
a very powerful tool to model high resolution observations 
like the ones from the TAURUS instrument (Walton et al. 
1990). 

EDITOR: PLACE FIGURE | HERE. 



4. Concluding remarks 

This paper illustrates the importance of 3D photoioniza- 
tion simulations when analyzing the various observational 
data currently obtained for PNe. Regarding the stellar tem- 
perature and luminosity, two different PNe are simulated. 
For each one, two density distributions are assumed, in or- 
der to test the constraints on the nebula geometry provided 
by the data. For each kind of observation (emission line 
spectroscopy, emission line imaging, emission line profiles) 
the results obtained with the two different geometries are 
compared. 

Regarding line emission from the whole nebula, the emis- 
sion line ratios do not discriminate between the two possible 
density distributions since the two models have the same 
average density, giving similar [S II] line intensity ratios. 
Emission line imaging showing an ellipsoidal shape with an 
intensity enhancement on the minor axis is obtained with 
both density distributions. Thus, this kind of observational 
data do not provide a method to distinguish between the 
two distributions. Density maps, obtained from emission 
line ratios that are density indicators, can provide one of the 



keys to the geometry puzzle, as shown in §3.4. The other 
key, related to the orientation of the axis of symmetry of 
the nebula, comes from the emission line profiles observed 
in different locations of the nebula. An asymmetric double- 
peaked profile is usually observed in PNe, but until now 
never theoretically reproduced by photoionization models. 



It can be generated assuming that the expansion velocity in- 
creases outward and that the axis of symmetry of the nebula 
makes an angle with the plane of the sky (see §4). The axis 
of symmetry tilt can also be visualized on the echelograms 
(§5). Notice that double-peaked profiles are also generated 
by pure hydrodynamic models (Mellema 1994). 

In brief, the different kinds of observations that are now 
available for planetary nebulae can only provide a complete 
insight of the object and its properties, if they are con- 
fronted with the results of photoionization simulations from 
a 3D code. 
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Fig. 1. — Emission- line imaging from models Al (top four 
panels) and A2 (botton four panels). For each model, H/3, 
[O III]A5007, He IIA4686, and [N II]A6583 intensity maps 
are shown. Both axes are in pixels. 

Fig. 2. — Same as Fig. [I] for models Bl (top four panels) 
and B2 (botton four panels). Both axes are in pixels. 

Fig. 3.— [S II] density map (from 200 to 900 cm" 3 ) for 
models Bl (top) and B2 (botton). The axes are in pixels. 

Fig. 4. — Size and positions of the aperture (top) used to 
obtain the [O III] velocity profiles shown in the botton nine 
panels for model Al; The nebula is rotated by 50° around 
the z axis and projected on the y direction. The axes are in 
pixels. The emission-line profiles (botton) are normalized. 

Fig. 5. — [O III] line image (top left), P-V diagrams ob- 
tained through centered horizontal and vertical slits (top 
right and botton left, respectively), and iso- velocity con- 
tours for V — 9 km/s. z and y are given in pixels. 



This piano tables was prepared with the AAS iATf^X macros v4.0. 
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Model 


Al 


A2 


Bl B2 


T e // (K) 


150000 


150000 


90000 90000 


L/L Q 


9260 


9260 


150 150 


njy(cav) (cm -3 ) 








400 400 


n H (shell) (cm -3 ) 


2000-4000 


3000 


550 - 1100 825 


R™ner/10 17 Cm 


1. 


0.5-1.5 


1. 0.5-1.5 


H/3 (10 35 erg/s) 


1.89 


1.88 


.048 .047 


HeII4686 


.330 


.328 


.065 .065 


[01] 6300+6363 


.227 


.232 


.402 .392 


[Oil] 3727 


1.83 


1.97 


2.72 2.63 


[OIII] 5007+4959 


24.0 


23.6 


6.72 6.92 


[OIII]4363 


.202 


.197 


.018 .019 


[Nil] 6583+6548 


1.88 


1.90 


3.08 2.99 


[Nil] 5755 


.032 


.032 


.034 .033 


[SIIJ6718 


.201 


.213 


.608 .598 


[SIIJ6732 


.297 


.307 


.624 .607 


T[o//7] (K) 


12173 


12122 


8558 8567 


< T e > 50 07 (K) 


12021 


11978 


8591 8596 


T[2vz/] (K) 


11621 


11565 


9568 9586 


< T e > 6583 (K) 


11194 


11225 


9582 9612 


n[s/j] (cm -3 ) 


2313 


2083 


657 629 


< n e > 67 i8 (cm -3 ) 


2491 


2272 


655 620 



Table 1: Input parameters: effective temperature and lu- 
minosity of the central ionizing black body, cavity density 
shell density, cavity radius or the semi-axis of the elliptical 
cavity (top). Emission line intensities relative to Hf3 (mid- 
dle). Physical quantities determined from line ratios and 



from average calculations, see § |2.5| (bottom) 
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Fig. 1.- 



H I 4861 [0 III] 5007 H I 4861 [0 III] 5007 
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Fig. 2.- 




Fig. 3.- 
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